Load packages
# rm(list = ls())
library(tictoc) # to monitor time
library(raster) # read raster files
library(rgdal) # Use GDAL functions
library(exactextractr) # Zonal statistics from raster
library(sf) # Working with shapefile layers
library(tidyverse) # data wrangling
library(plyr) # data wrangling
library(dplyr) # data wrangling
library(ggpmisc) # To place geom_text within plot
library(plotly) # Interactive visualization
library(patchwork) # This package for arranging multiple ggplots
print("Packages loaded successfully!")
## [1] "Packages loaded successfully!"
Functions specific to this RMD
#########################################
# Creating function for exponential fit
myexpfit <- function(response_var, pred_var, mydata){
# print(paste0("Model = ", "log(",response_var, ")~", pred_var))
model.linear <- lm(paste0("log(",response_var, ")~", pred_var), data = mydata)
a1 <- exp(model.linear$coefficients[1])
b1 <- (model.linear$coefficients[2])
R2 <- summary(model.linear)$r.squared
myresults <- c(a1, b1, R2)
return(myresults)
}
#########################################
# My GGPLOT for yield map
my_yield_plot <- function(yield_df){
fieldtitle <- levels(yield_df$Fieldname)
# Legend labeling - Splitting the yield range into equal classes
minyld <- round(min(yield_df$Yield_Mg.ha),0) # minimum yield
maxyld <- round(max(yield_df$Yield_Mg.ha),0) # maximum yield
Nclasses <- 4 # Number of yield classes to split from the whole yield range
mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
YE_colors <- c("blue","cyan","lightgreen","yellow","orange", "red")
yld_plot <- ggplot(yield_df, aes(x = Longitude, y = Latitude)) +
geom_point(aes(color = Yield_Mg.ha), size = 0.7) +
theme_bw() +
# facet_grid(Field~Year, scales = "free") +
theme(aspect.ratio = 1.0) +
# ggtitle("Actual yield (Mg/ha)")+
facet_wrap(~Fieldname, scales = "free") +
# theme(legend.position = c(0.85, 0.27))+
scale_colour_gradientn(colours = YE_colors,
name = "Yield (Mg/ha)",
limits=c(minyld,maxyld),
breaks = mybreaks,
labels = round(mybreaks, 0)) +
guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
yld_plot
}
#########################################
my_pred_plot <- function(sub_dat){
fieldtitle <- levels(sub_dat$Fieldname)
# Legend labeling - Splitting the yield range into equal classes
minyld <- round(min(sub_dat$Pred_Yield),0) # minimum yield
maxyld <- round(max(sub_dat$Pred_Yield),0) # maximum yield
Nclasses <- 4 # Number of yield classes to split from the whole yield range
mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
YE_colors <- c("blue","cyan","lightgreen","yellow","orange", "red")
pred_map <- ggplot(sub_dat, aes(x = Longitude, y = Latitude)) +
geom_point(aes(color = Pred_Yield), size = 0.7) +
theme_bw() +
# facet_grid(Field~Year, scales = "free") +
theme(aspect.ratio = 1.0) +
# ggtitle("Predicted yield (Mg/ha)")+
facet_wrap(~Fieldname, scales = "free") +
# theme(legend.position = c(0.85, 0.27))+
scale_colour_gradientn(colours = YE_colors,
name = "Yield (Mg/ha)",
limits=c(minyld,maxyld),
breaks = mybreaks,
labels = round(mybreaks, 0)) +
guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
pred_map
}
#########################################
## Actual and predicted yield map
make_yield_maps <- function(sub_datf){
sel_dat4 <- sub_datf %>%
select(Fieldname, Latitude, Longitude, Yield_Mg.ha, Pred_Yield) %>%
gather("yield_type", "Value", 4:5) %>%
droplevels()
sel_dat4$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
sel_dat4$yield_type <- factor(sel_dat4$yield_type, levels = c("Actual yield", "Predicted yield"))
minyld <- round(min(sel_dat4$Value),0) # minimum yield
maxyld <- round(max(sel_dat4$Value),0) # maximum yield
Nclasses <- 4 # Number of yield classes to split from the whole yield range
mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
YE_colors <- c("blue","cyan","lightgreen","yellow","orange", "red")
yld_plot <- ggplot(sel_dat4, aes(x = Longitude, y = Latitude)) +
geom_point(aes(color = Value), size = 0.7) +
theme_bw() +
# facet_grid(Field~Year, scales = "free") +
theme(aspect.ratio = 1.0) +
# ggtitle("Actual yield (Mg/ha)")+
facet_wrap(~yield_type, scales = "free") +
# theme(legend.position = c(0.85, 0.27))+
scale_colour_gradientn(colours = YE_colors,
name = "Yield (Mg/ha)",
limits=c(minyld,maxyld),
breaks = mybreaks,
labels = round(mybreaks, 0)) +
guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
yld_plot
}
Load dataset
load("Quantix_deep_learning.RData", verbose = TRUE)
## Loading objects:
## dat
str(dat)
## 'data.frame': 1270534 obs. of 46 variables:
## $ Fieldname : Factor w/ 8 levels "DMD_GH1","PSF_111",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ FlightDate : Factor w/ 22 levels "2019/06/07","2019/06/09",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ DOY : int 158 158 158 158 158 158 158 158 158 158 ...
## $ Week : Factor w/ 14 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Latitude : num 43 43 43 43 43 ...
## $ Longitude : num -76.6 -76.6 -76.6 -76.6 -76.6 ...
## $ Yield : num 182 188 191 194 197 ...
## $ Yield_Mg.ha : num 11.4 11.8 12 12.2 12.4 ...
## $ Type : Factor w/ 2 levels "grain","silage": 1 1 1 1 1 1 1 1 1 1 ...
## $ N : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ Strip : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rgn_red : num 18 18 18 18 18 18 18 18 18 18 ...
## $ rgn_green : num 16 15 15 15 15 15 15 15 15 16 ...
## $ rgn_nir : num 36 36 35 35 35 35 35 35 35 35 ...
## $ rgb_red : num 115 115 116 117 118 118 119 118 118 119 ...
## $ rgb_green : num 89 89 90 90 91 91 92 91 91 92 ...
## $ rgb_blue : num 67 67 68 68 69 69 70 69 69 70 ...
## $ NDVI : num 0.333 0.333 0.321 0.321 0.321 ...
## $ GNDVI : num 0.385 0.412 0.4 0.4 0.4 ...
## $ EVI2 : num 0.597 0.616 0.59 0.59 0.59 ...
## $ SR : num 2 2 1.94 1.94 1.94 ...
## $ EXG : num -4 -4 -4 -5 -5 -5 -5 -5 -5 -5 ...
## $ TGI : num 3.28 3.28 3.28 2.89 2.89 2.89 2.89 2.89 2.89 2.89 ...
## $ Mgmt_zone : Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 3 2 2 2 2 ...
## $ Climod_GDD : int 10 10 10 10 10 10 10 10 10 10 ...
## $ Climod_Temp_C : num 15.6 15.6 15.6 15.6 15.6 ...
## $ Climod_Precip_mm : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ACIS_GDD : int 10 10 10 10 10 10 10 10 10 10 ...
## $ ACIS_Temp_C : num 15.6 15.6 15.6 15.6 15.6 ...
## $ ACIS_Precip_mm : num 0 0 0 0 0 0 0 0 0 0 ...
## $ NASA_PAR : num 159 159 159 159 159 ...
## $ NASA_Temp_C : num 17.3 17.3 17.3 17.3 17.3 ...
## $ NASA_SH : num 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 7.2 ...
## $ NASA_RH : num 59.7 59.7 59.7 59.7 59.7 ...
## $ NASA_Precip_mm : num 0 0 0 0 0 0 0 0 0 0 ...
## $ NASA_WS : num 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 0.56 ...
## $ NASA_SurfWet : num 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
## $ NASA_RootWet : num 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 ...
## $ NASA_ProfWet : num 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 0.66 ...
## $ DEM : int 125 125 125 125 125 125 125 125 125 125 ...
## $ DEM_type_Flat_Area : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ DEM_type_Lower_Slope : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 2 2 ...
## $ DEM_type_Middle_Slope: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ DEM_type_Ridge : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 1 1 1 ...
## $ DEM_type_Upper_Slope : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ DEM_type_Valley : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
r2df_all <- data.frame()
Data preprocessing
Check for missing values
# Check if any columns have missing values
sapply(dat, anyNA)
## Fieldname FlightDate DOY
## FALSE FALSE FALSE
## Week Latitude Longitude
## FALSE FALSE FALSE
## Yield Yield_Mg.ha Type
## FALSE FALSE FALSE
## N Strip rgn_red
## TRUE TRUE FALSE
## rgn_green rgn_nir rgb_red
## FALSE FALSE FALSE
## rgb_green rgb_blue NDVI
## FALSE FALSE FALSE
## GNDVI EVI2 SR
## FALSE FALSE FALSE
## EXG TGI Mgmt_zone
## FALSE FALSE TRUE
## Climod_GDD Climod_Temp_C Climod_Precip_mm
## FALSE FALSE FALSE
## ACIS_GDD ACIS_Temp_C ACIS_Precip_mm
## FALSE FALSE FALSE
## NASA_PAR NASA_Temp_C NASA_SH
## FALSE FALSE FALSE
## NASA_RH NASA_Precip_mm NASA_WS
## FALSE FALSE FALSE
## NASA_SurfWet NASA_RootWet NASA_ProfWet
## FALSE FALSE FALSE
## DEM DEM_type_Flat_Area DEM_type_Lower_Slope
## FALSE FALSE FALSE
## DEM_type_Middle_Slope DEM_type_Ridge DEM_type_Upper_Slope
## FALSE FALSE FALSE
## DEM_type_Valley
## FALSE
# Names of columns that has missing values
names(dat)[sapply(dat, anyNA)]
## [1] "N" "Strip" "Mgmt_zone"
Plot all yield data
g_dat <- dat %>%
dplyr::filter(Type == "grain") %>%
droplevels()
my_yield_plot(g_dat)
s_dat <- dat %>%
dplyr::filter(Type == "silage") %>%
droplevels()
my_yield_plot(s_dat)
Exploratory data analysis
Distribution plots
# Violin plot
vplt <- ggplot(dat, aes(x = Fieldname, y = Yield_Mg.ha, fill = Fieldname)) +
geom_violin(alpha = 0.4) +
theme_bw() +
facet_wrap(~Type, nrow = 2, scales = "free_y")
vplt
dplt <- ggplot(dat, aes(x = Yield_Mg.ha, y = ..scaled.., fill = Fieldname)) +
geom_density(alpha = 0.4, size = 1) +
theme_bw() +
facet_wrap(~Type, nrow = 2, scales = "free")
dplt
GPS points vs yield
# "DMD_GH1" "PSF_111" "PSF_12" "SLS_ABH" "SLS_NS" "SSF_121" "SSF_202" "SSF_66"
# selfield <- "DMD_GH1"
sub_datf <- dat %>%
filter(Week == 1 & Type == "grain") %>%
droplevels()
# Legend labeling - Splitting the yield range into equal classes
minyld <- round(min(sub_datf$Yield_Mg.ha),0) # minimum yield
maxyld <- round(max(sub_datf$Yield_Mg.ha),0) # maximum yield
Nclasses <- 4 # Number of yield classes to split from the whole yield range
mybreaks <- round(seq(minyld, maxyld, len = (Nclasses+1)),0)
YE_colors <- c("blue","cyan","lightgreen","yellow","orange", "red")
#=====================
lat_plot <- ggplot(sub_datf, aes(x = Latitude, y = Yield_Mg.ha)) +
geom_point(aes(color = Yield_Mg.ha), size = 0.7) +
theme_bw() +
# facet_grid(Field~Year, scales = "free") +
theme(aspect.ratio = 1.0) +
ggtitle("Latitude vs Yield")+
facet_wrap(~Fieldname, scales = "free") +
# theme(legend.position = c(0.85, 0.27))+
scale_colour_gradientn(colours = YE_colors,
name = "Yield (Mg/ha)",
limits=c(minyld,maxyld),
breaks = mybreaks,
labels = round(mybreaks, 0)) +
guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
lat_plot
#=====================
long_plot <- ggplot(sub_datf, aes(x = Longitude, y = Yield_Mg.ha)) +
geom_point(aes(color = Yield_Mg.ha), size = 0.7) +
theme_bw() +
# facet_grid(Field~Year, scales = "free") +
theme(aspect.ratio = 1.0) +
ggtitle("Longitude vs Yield")+
facet_wrap(~Fieldname, scales = "free") +
# theme(legend.position = c(0.85, 0.27))+
scale_colour_gradientn(colours = YE_colors,
name = "Yield (Mg/ha)",
limits=c(minyld,maxyld),
breaks = mybreaks,
labels = round(mybreaks, 0)) +
guides(colour = guide_colourbar(barheight = unit(1.8, "in")))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
long_plot
Approach 1
1.a. Use only strip data for exponential model
# "DMD_GH1" "PSF_111" "PSF_12" "SLS_ABH" "SLS_NS" "SSF_121" "SSF_202" "SSF_66"
# "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"
# Best grain field - DMD_GH1
# Best silage filed - SSF_121, PSF_12
vi <- c("NDVI", "GNDVI", "EVI2", "SR", "EXG", "TGI")
selfield <- "DMD_GH1"
# seltype <- "silage"
sub_dat <- dat %>%
filter(!is.na(N) & Fieldname == selfield) %>%
droplevels()
my_yield_plot(sub_dat)
# "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"
lvl <- levels(sub_dat$Week)
sub_dat$Strip <- as.numeric(sub_dat$Strip)
#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(sub_dat, aes(x = NDVI, y = Yield_Mg.ha)) +
geom_point(aes(color = Week)) +
theme_bw() +
theme(aspect.ratio = 1) +
ggtitle("NDVI vs Yield") +
xlab("NDVI") +
ylab("Yield (Mg/ha)") +
xlim(c(0.0, 1.0))
sct_plt
#==============
r2df <- data.frame()
for (i in 1:length(lvl)){
# print(i)
wsub_dat <- sub_dat %>%
filter(Week == lvl[i]) %>%
droplevels()
for (j in 1:length(vi)){
dep_var <- "Yield_Mg.ha"
ind_var <- vi[j]
# cat(dep_var, "---", ind_var, "\n")
mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
## RMSE calculation ======
a_fit <- as.numeric(mod[1])
b_fit <- as.numeric(mod[2])
rsq <- as.numeric(mod[3])
ssub_dat <- dat %>%
filter(Week == lvl[i] & Fieldname == selfield) %>%
droplevels()
ssub_dat$Pred_Yield <- a_fit*exp(b_fit*ssub_dat[[ind_var]])
ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
sqe <- (ssub_dat$error^2)
mse <- mean(sqe)
rmse_week <- sqrt(mse)
# cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
# cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
## ========================
r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
r2df <- rbind(r2df, r2)
}
# print(paste0(i, " = ", mod[3]))
}
colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")
num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame': 78 obs. of 6 variables:
## $ Weeks: num 1 1 1 1 1 1 2 2 2 2 ...
## $ a : num 11.13 10.62 11.21 9.89 14.35 ...
## $ b : num 0.87957 0.84782 0.47224 0.20573 -0.00694 ...
## $ R2 : num 0.0329 0.0314 0.0371 0.0349 0.0249 ...
## $ RMSE : num 2.29 2.21 2.3 3.85 2.02 ...
## $ VI : chr "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("1a", length(lvl))
r2df_all <- rbind(r2df_all, r2df)
1.b. Use lower and upper 10% of strip data
# "DMD_GH1" "PSF_111" "PSF_12" "SLS_ABH" "SLS_NS" "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"
sub_dat <- dat %>%
filter(!is.na(N) & Fieldname == selfield) %>%
droplevels()
my_yield_plot(sub_dat)
# Evaluate qunatile cutoffs
pert <- 0.10
cutoffs <- quantile(sub_dat$Yield_Mg.ha, c(pert, (1-pert)))
cat("Cutoff for field", selfield, " =", cutoffs)
## Cutoff for field SSF_66 = 59.6752 81.1395
# Subset points within cutoff
qsub_dat <- sub_dat %>%
filter(Yield_Mg.ha <= cutoffs[[1]] | Yield_Mg.ha >= cutoffs[[2]]) %>%
droplevels()
my_yield_plot(qsub_dat)
#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(qsub_dat, aes(x = NDVI, y = Yield_Mg.ha)) +
geom_point(aes(color = Week)) +
theme_bw() +
theme(aspect.ratio = 1) +
xlab("NDVI") +
ylab("Yield (Mg/ha)") +
xlim(c(0.0, 1.0))
sct_plt
#==============
lvl <- levels(qsub_dat$Week)
r2df <- data.frame()
for (i in 1:length(lvl)){
# print(i)
wsub_dat <- qsub_dat %>%
filter(Week == lvl[i]) %>%
droplevels()
for (j in 1:length(vi)){
dep_var <- "Yield_Mg.ha"
ind_var <- vi[j]
# cat(dep_var, "---", ind_var, "\n")
mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
## RMSE calculation ======
a_fit <- as.numeric(mod[1])
b_fit <- as.numeric(mod[2])
rsq <- as.numeric(mod[3])
ssub_dat <- dat %>%
filter(Week == lvl[i] & Fieldname == selfield) %>%
droplevels()
ssub_dat$Pred_Yield <- a_fit*exp(b_fit*ssub_dat[[ind_var]])
ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
sqe <- (ssub_dat$error^2)
mse <- mean(sqe)
rmse_week <- sqrt(mse)
# cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
# cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
## ========================
r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
r2df <- rbind(r2df, r2)
}
# print(paste0(i, " = ", mod[3]))
}
colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")
num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame': 78 obs. of 6 variables:
## $ Weeks: num 1 1 1 1 1 1 2 2 2 2 ...
## $ a : num 82.6 101.3 83.4 80.3 67.1 ...
## $ b : num -0.64658 -1.11741 -0.38113 -0.08987 0.00485 ...
## $ R2 : num 0.0106 0.0228 0.01513 0.01428 0.00636 ...
## $ RMSE : num 8.97 9 8.97 8.97 9.16 ...
## $ VI : chr "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("1b", length(lvl))
r2df_all <- rbind(r2df_all, r2df)
Approach 2
2.a. Use whole field data without strip –> is.na(N)
# "DMD_GH1" "PSF_111" "PSF_12" "SLS_ABH" "SLS_NS" "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"
sub_dat <- dat %>%
filter(is.na(N) & Fieldname == selfield) %>%
droplevels()
my_yield_plot(sub_dat)
#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(sub_dat, aes(x = NDVI, y = Yield_Mg.ha)) +
geom_point(aes(color = Week)) +
theme_bw() +
theme(aspect.ratio = 1) +
xlab("NDVI") +
ylab("Yield (Mg/ha)") +
xlim(c(0.0, 1.0))
sct_plt
#==============
lvl <- levels(sub_dat$Week)
r2df <- data.frame()
for (i in 1:length(lvl)){
# print(i)
wsub_dat <- sub_dat %>%
filter(Week == lvl[i]) %>%
droplevels()
for (j in 1:length(vi)){
dep_var <- "Yield_Mg.ha"
ind_var <- vi[j]
# cat(dep_var, "---", ind_var, "\n")
mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
## RMSE calculation ======
a_fit <- as.numeric(mod[1])
b_fit <- as.numeric(mod[2])
rsq <- as.numeric(mod[3])
ssub_dat <- dat %>%
filter(Week == lvl[i] & Fieldname == selfield) %>%
droplevels()
ssub_dat$Pred_Yield <- a_fit*exp(b_fit*ssub_dat[[ind_var]])
ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
sqe <- (ssub_dat$error^2)
mse <- mean(sqe)
rmse_week <- sqrt(mse)
# cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
# cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
## ========================
r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
r2df <- rbind(r2df, r2)
}
# print(paste0(i, " = ", mod[3]))
}
colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")
num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame': 78 obs. of 6 variables:
## $ Weeks: num 1 1 1 1 1 1 2 2 2 2 ...
## $ a : num 94.2 101 90.9 85.6 68.5 ...
## $ b : num -1.06597 -1.09442 -0.5315 -0.12254 -0.00442 ...
## $ R2 : num 0.0535 0.0397 0.0523 0.0393 0.0132 ...
## $ RMSE : num 8.96 8.96 8.95 8.94 9.09 ...
## $ VI : chr "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("2a", length(lvl))
r2df_all <- rbind(r2df_all, r2df)
2.b. Use lower and upper 10% of whole field data
# "DMD_GH1" "PSF_111" "PSF_12" "SLS_ABH" "SLS_NS" "SSF_121" "SSF_202" "SSF_66"
selfield <- "SSF_66"
sub_dat <- dat %>%
filter(is.na(N) & Fieldname == selfield) %>%
droplevels()
my_yield_plot(sub_dat)
# Evaluate quantile cutoffs
pert <- 0.10
cutoffs <- quantile(sub_dat$Yield_Mg.ha, c(pert, (1-pert)))
cat("Cutoff for field", selfield, " =", cutoffs)
## Cutoff for field SSF_66 = 56.8841 81.4359
# Subset points within cutoff
qsub_dat <- sub_dat %>%
filter(Yield_Mg.ha <= cutoffs[[1]] | Yield_Mg.ha >= cutoffs[[2]]) %>%
droplevels()
my_yield_plot(qsub_dat)
#==============
# Scatter plot between NDVI and Yield
sct_plt <- ggplot(qsub_dat, aes(x = NDVI, y = Yield_Mg.ha)) +
geom_point(aes(color = Week)) +
theme_bw() +
theme(aspect.ratio = 1) +
xlab("NDVI") +
ylab("Yield (Mg/ha)") +
xlim(c(0.0, 1.0))
sct_plt
#==============
lvl <- levels(qsub_dat$Week)
r2df <- data.frame()
for (i in 1:length(lvl)){
# print(i)
wsub_dat <- qsub_dat %>%
filter(Week == lvl[i]) %>%
droplevels()
for (j in 1:length(vi)){
dep_var <- "Yield_Mg.ha"
ind_var <- vi[j]
# cat(dep_var, "---", ind_var, "\n")
mod <- myexpfit("Yield_Mg.ha", ind_var, wsub_dat)
## RMSE calculation ======
a_fit <- as.numeric(mod[1])
b_fit <- as.numeric(mod[2])
rsq <- as.numeric(mod[3])
ssub_dat <- dat %>%
filter(Week == lvl[i] & Fieldname == selfield) %>%
droplevels()
ssub_dat$Pred_Yield <- a_fit*exp(b_fit*ssub_dat[[ind_var]])
ssub_dat$error <- (ssub_dat$Pred_Yield - ssub_dat$Yield_Mg.ha)
sqe <- (ssub_dat$error^2)
mse <- mean(sqe)
rmse_week <- sqrt(mse)
cat(lvl[i], ") a_fit = ", a_fit, "; b_fit = ", b_fit, ";\t R2 = ", rsq,";\t RMSE = ", rmse_week, "\n")
# cat("RMSE for", selfield, "on week", lvl[i], "is", rmse, "Mg/ha")
## ========================
r2 <- c(as.numeric(c(lvl[i], mod[1], mod[2], mod[3], rmse_week)), ind_var)
r2df <- rbind(r2df, r2)
}
# print(paste0(i, " = ", mod[3]))
}
## 1 ) a_fit = 114.5396 ; b_fit = -1.775352 ; R2 = 0.1047697 ; RMSE = 9.3737
## 1 ) a_fit = 127.595 ; b_fit = -1.814529 ; R2 = 0.07242787 ; RMSE = 9.418363
## 1 ) a_fit = 104.0244 ; b_fit = -0.8232189 ; R2 = 0.0940588 ; RMSE = 9.351495
## 1 ) a_fit = 85.56061 ; b_fit = -0.138763 ; R2 = 0.05338813 ; RMSE = 9.445993
## 1 ) a_fit = 69.44581 ; b_fit = -0.01730106 ; R2 = 0.08496018 ; RMSE = 9.968001
## 1 ) a_fit = 80.12686 ; b_fit = -0.02549843 ; R2 = 0.0559519 ; RMSE = 9.91363
## 2 ) a_fit = 97.68063 ; b_fit = -0.8558484 ; R2 = 0.08915737 ; RMSE = 9.710312
## 2 ) a_fit = 100.4366 ; b_fit = -0.8612361 ; R2 = 0.05955035 ; RMSE = 9.634497
## 2 ) a_fit = 90.54516 ; b_fit = -0.3716412 ; R2 = 0.08503231 ; RMSE = 9.616456
## 2 ) a_fit = 75.43675 ; b_fit = -0.04808592 ; R2 = 0.08371296 ; RMSE = 9.359314
## 2 ) a_fit = 72.97827 ; b_fit = -0.01113459 ; R2 = 0.07962387 ; RMSE = 10.06643
## 2 ) a_fit = 77.92035 ; b_fit = -0.01842857 ; R2 = 0.05551291 ; RMSE = 9.937669
## 3 ) a_fit = 130.2284 ; b_fit = -1.772073 ; R2 = 0.08667496 ; RMSE = 9.545287
## 3 ) a_fit = 151.7014 ; b_fit = -1.912627 ; R2 = 0.07356775 ; RMSE = 9.665523
## 3 ) a_fit = 120.4115 ; b_fit = -0.8546063 ; R2 = 0.08906408 ; RMSE = 9.555316
## 3 ) a_fit = 93.15247 ; b_fit = -0.1546254 ; R2 = 0.05405708 ; RMSE = 9.430931
## 3 ) a_fit = 77.07453 ; b_fit = -0.02416504 ; R2 = 0.07798439 ; RMSE = 9.901467
## 3 ) a_fit = 87.24738 ; b_fit = -0.03597171 ; R2 = 0.05820835 ; RMSE = 10.1733
## 4 ) a_fit = 39.96259 ; b_fit = 0.7992891 ; R2 = 0.049653 ; RMSE = 9.760227
## 4 ) a_fit = 35.31187 ; b_fit = 0.9765932 ; R2 = 0.05753524 ; RMSE = 9.731614
## 4 ) a_fit = 44.29221 ; b_fit = 0.3180522 ; R2 = 0.04569029 ; RMSE = 9.737102
## 4 ) a_fit = 55.23613 ; b_fit = 0.03822832 ; R2 = 0.0205722 ; RMSE = 9.777046
## 4 ) a_fit = 80.0245 ; b_fit = -0.01120671 ; R2 = 0.02603877 ; RMSE = 9.936435
## 4 ) a_fit = 95.7974 ; b_fit = -0.03033375 ; R2 = 0.05545126 ; RMSE = 9.916124
## 5 ) a_fit = 21.26627 ; b_fit = 1.377981 ; R2 = 0.1782023 ; RMSE = 9.033781
## 5 ) a_fit = 20.33525 ; b_fit = 1.457281 ; R2 = 0.1927672 ; RMSE = 9.040896
## 5 ) a_fit = 26.53846 ; b_fit = 0.5126518 ; R2 = 0.2002983 ; RMSE = 9.111857
## 5 ) a_fit = 42.33944 ; b_fit = 0.04008052 ; R2 = 0.1821739 ; RMSE = 9.767657
## 5 ) a_fit = 50.62741 ; b_fit = 0.007041083 ; R2 = 0.05393396 ; RMSE = 9.303907
## 5 ) a_fit = 52.54574 ; b_fit = 0.01057815 ; R2 = 0.04159836 ; RMSE = 9.362912
## 6 ) a_fit = 17.90351 ; b_fit = 1.559624 ; R2 = 0.153909 ; RMSE = 9.08576
## 6 ) a_fit = 15.60691 ; b_fit = 1.747906 ; R2 = 0.163055 ; RMSE = 9.040844
## 6 ) a_fit = 23.14956 ; b_fit = 0.5700638 ; R2 = 0.1629603 ; RMSE = 9.067256
## 6 ) a_fit = 40.48376 ; b_fit = 0.04118642 ; R2 = 0.1547518 ; RMSE = 9.504802
## 6 ) a_fit = 75.83658 ; b_fit = -0.004390261 ; R2 = 0.009154759 ; RMSE = 9.77616
## 6 ) a_fit = 84.54597 ; b_fit = -0.01349674 ; R2 = 0.02717307 ; RMSE = 9.746971
## 7 ) a_fit = 12.22936 ; b_fit = 1.892968 ; R2 = 0.1274779 ; RMSE = 9.175553
## 7 ) a_fit = 10.66996 ; b_fit = 2.097191 ; R2 = 0.1527872 ; RMSE = 9.029418
## 7 ) a_fit = 18.26674 ; b_fit = 0.6438008 ; R2 = 0.1446185 ; RMSE = 9.052002
## 7 ) a_fit = 43.57484 ; b_fit = 0.0227518 ; R2 = 0.1088754 ; RMSE = 9.814816
## 7 ) a_fit = 84.75155 ; b_fit = -0.007686897 ; R2 = 0.03953933 ; RMSE = 10.2967
## 7 ) a_fit = 85.4768 ; b_fit = -0.01445088 ; R2 = 0.04640809 ; RMSE = 10.3253
## 8 ) a_fit = 7.357011 ; b_fit = 2.434166 ; R2 = 0.1423155 ; RMSE = 8.993199
## 8 ) a_fit = 7.491671 ; b_fit = 2.466606 ; R2 = 0.1648986 ; RMSE = 8.95643
## 8 ) a_fit = 14.61957 ; b_fit = 0.7389926 ; R2 = 0.1518354 ; RMSE = 8.95934
## 8 ) a_fit = 36.13157 ; b_fit = 0.03005452 ; R2 = 0.1783001 ; RMSE = 9.1046
## 8 ) a_fit = 109.8534 ; b_fit = -0.01710957 ; R2 = 0.2021991 ; RMSE = 9.39526
## 8 ) a_fit = 108.4526 ; b_fit = -0.03030211 ; R2 = 0.2090256 ; RMSE = 9.358248
## 9 ) a_fit = 8.986937 ; b_fit = 2.206359 ; R2 = 0.1272976 ; RMSE = 9.070061
## 9 ) a_fit = 7.650902 ; b_fit = 2.440692 ; R2 = 0.1586471 ; RMSE = 9.021347
## 9 ) a_fit = 14.87923 ; b_fit = 0.7285043 ; R2 = 0.1469295 ; RMSE = 9.019683
## 9 ) a_fit = 37.4024 ; b_fit = 0.02766957 ; R2 = 0.1659178 ; RMSE = 9.36852
## 9 ) a_fit = 99.99953 ; b_fit = -0.01418411 ; R2 = 0.1435705 ; RMSE = 9.552597
## 9 ) a_fit = 100.1835 ; b_fit = -0.02569391 ; R2 = 0.1505604 ; RMSE = 9.516102
## 10 ) a_fit = 11.57471 ; b_fit = 1.960843 ; R2 = 0.1132833 ; RMSE = 9.217313
## 10 ) a_fit = 10.2941 ; b_fit = 2.127497 ; R2 = 0.1355196 ; RMSE = 9.063766
## 10 ) a_fit = 17.07697 ; b_fit = 0.6752103 ; R2 = 0.1347896 ; RMSE = 9.052663
## 10 ) a_fit = 43.64423 ; b_fit = 0.02357225 ; R2 = 0.09105606 ; RMSE = 9.576446
## 10 ) a_fit = 80.67912 ; b_fit = -0.006875877 ; R2 = 0.05818177 ; RMSE = 9.600197
## 10 ) a_fit = 80.468 ; b_fit = -0.01210274 ; R2 = 0.05960173 ; RMSE = 9.592897
## 11 ) a_fit = 3.369919 ; b_fit = 3.344318 ; R2 = 0.1408751 ; RMSE = 9.132297
## 11 ) a_fit = 3.217137 ; b_fit = 3.458154 ; R2 = 0.2017893 ; RMSE = 9.003591
## 11 ) a_fit = 8.643776 ; b_fit = 1.013074 ; R2 = 0.1740437 ; RMSE = 9.007777
## 11 ) a_fit = 41.28896 ; b_fit = 0.02641921 ; R2 = 0.1010583 ; RMSE = 9.485659
## 11 ) a_fit = 105.4118 ; b_fit = -0.01461405 ; R2 = 0.1917603 ; RMSE = 9.27535
## 11 ) a_fit = 103.7901 ; b_fit = -0.02526568 ; R2 = 0.1948443 ; RMSE = 9.245671
## 12 ) a_fit = 2.8717 ; b_fit = 3.513259 ; R2 = 0.153327 ; RMSE = 9.169528
## 12 ) a_fit = 4.033029 ; b_fit = 3.199581 ; R2 = 0.1682395 ; RMSE = 9.06045
## 12 ) a_fit = 9.612162 ; b_fit = 0.9600938 ; R2 = 0.154464 ; RMSE = 9.040218
## 12 ) a_fit = 40.85059 ; b_fit = 0.02620225 ; R2 = 0.1183588 ; RMSE = 9.714413
## 12 ) a_fit = 92.04569 ; b_fit = -0.01179612 ; R2 = 0.1512497 ; RMSE = 9.75852
## 12 ) a_fit = 91.31656 ; b_fit = -0.02042833 ; R2 = 0.1523465 ; RMSE = 9.732931
## 13 ) a_fit = 13.97952 ; b_fit = 1.765124 ; R2 = 0.113335 ; RMSE = 9.107881
## 13 ) a_fit = 12.26431 ; b_fit = 2.01955 ; R2 = 0.1525834 ; RMSE = 9.115751
## 13 ) a_fit = 18.36556 ; b_fit = 0.6738635 ; R2 = 0.1539626 ; RMSE = 9.050074
## 13 ) a_fit = 47.63929 ; b_fit = 0.01957556 ; R2 = 0.08607894 ; RMSE = 9.465863
## 13 ) a_fit = 77.01419 ; b_fit = -0.004729239 ; R2 = 0.04437182 ; RMSE = 9.686164
## 13 ) a_fit = 77.0385 ; b_fit = -0.008274721 ; R2 = 0.04717384 ; RMSE = 9.672191
colnames(r2df) <- c("Weeks", "a", "b", "R2", "RMSE", "VI")
num_cols <- c("Weeks", "a", "b", "R2", "RMSE")
r2df[num_cols] <- lapply(r2df[num_cols], as.numeric)
str(r2df)
## 'data.frame': 78 obs. of 6 variables:
## $ Weeks: num 1 1 1 1 1 1 2 2 2 2 ...
## $ a : num 114.5 127.6 104 85.6 69.4 ...
## $ b : num -1.7754 -1.8145 -0.8232 -0.1388 -0.0173 ...
## $ R2 : num 0.1048 0.0724 0.0941 0.0534 0.085 ...
## $ RMSE : num 9.37 9.42 9.35 9.45 9.97 ...
## $ VI : chr "NDVI" "GNDVI" "EVI2" "SR" ...
r2df$Fieldname <- rep(selfield, length(lvl))
r2df$Approach <- rep("2b", length(lvl))
r2df_all <- rbind(r2df_all, r2df)
All vegetation indices + All weeks + All fields
load("Exponential_results_March31.RData", verbose = TRUE)
## Loading objects:
## results_df
Plotting performance data - All VI + All weeks + All fields + All approaches
# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"
# selfield <- "DMD_GH1"
# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
results_df$VI <- factor(results_df$VI, levels = c("TGI", "EXG", "SR", "GNDVI", "EVI2", "NDVI"))
# selVI <- "EVI2"
sub_results <- results_df %>%
# filter(Fieldname == selfield) %>%
droplevels()
bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
# geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7) +
geom_bar(stat="identity", position=position_dodge(), width = 0.6) +
theme_bw() +
theme(legend.position = "right", legend.box = "horizontal") +
# theme(legend.position = "none") +
# guides(colour = guide_legend(nrow = 1)) +
# ggtitle(selfield) +
# facet_wrap(~Approach) +
facet_grid(Fieldname~Approach) +
# facet_grid(Fieldname~Approach) +
scale_x_continuous(breaks = seq(1, 14, 1)) +
# scale_y_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
# geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2, color="black",
# position = position_dodge(0), size=2, angle = 90)+
# scale_fill_brewer(palette="Reds") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, face = "bold"),
# axis.title.y = element_blank(),
axis.text = element_text(size = 8, color = "black"),
# axis.text.y = element_blank(),
# axis.ticks.length=unit(-1.5, "mm"),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
bar_plot
# ggplotly(bar_plot)
# ggsave("All_VI_performance.pdf", width = 12, height = 8, units = "in")
Plotting performance data - NDVI
# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"
# selfield <- "SSF_66"
# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))
selVI <- "NDVI"
sub_results <- results_df %>%
filter(VI == selVI) %>%
droplevels()
bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
theme_bw() +
# theme(legend.position = "bottom", legend.box = "horizontal") +
theme(legend.position = "none") +
guides(colour = guide_legend(nrow = 1)) +
ggtitle(selVI) +
facet_grid(Fieldname~Approach) +
scale_x_continuous(breaks = seq(1, 14, 1)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2, color="black",
position = position_dodge(0), size=2, angle = 90)+
# scale_fill_brewer(palette="Reds") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, face = "bold"),
# axis.title.y = element_blank(),
axis.text = element_text(size = 8, color = "black"),
# axis.text.y = element_blank(),
# axis.ticks.length=unit(-1.5, "mm"),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
bar_plot
# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")
Plotting performance data - EVI2
# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"
# selfield <- "SSF_66"
# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))
selVI <- "EVI2"
sub_results <- results_df %>%
filter(VI == selVI) %>%
droplevels()
bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
theme_bw() +
# theme(legend.position = "bottom", legend.box = "horizontal") +
theme(legend.position = "none") +
guides(colour = guide_legend(nrow = 1)) +
ggtitle(selVI) +
facet_grid(Fieldname~Approach) +
scale_x_continuous(breaks = seq(1, 14, 1)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2, color="black",
position = position_dodge(0), size=2, angle = 90)+
# scale_fill_brewer(palette="Reds") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, face = "bold"),
# axis.title.y = element_blank(),
axis.text = element_text(size = 8, color = "black"),
# axis.text.y = element_blank(),
# axis.ticks.length=unit(-1.5, "mm"),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
bar_plot
# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")
Plotting performance data - GNDVI
# "DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121" "SSF_66", "SSF_202", "PSF_111", "PSF_12"
# selfield <- "SSF_66"
# results_df$yield_type <- recode_factor(factor(sel_dat4$yield_type), Pred_Yield = "Predicted yield", Yield_Mg.ha = "Actual yield")
results_df$Fieldname <- factor(results_df$Fieldname, levels = c("DMD_GH1", "SLS_ABH", "SLS_NS", "SSF_121", "SSF_66", "SSF_202", "PSF_111", "PSF_12"))
selVI <- "GNDVI"
sub_results <- results_df %>%
filter(VI == selVI) %>%
droplevels()
bar_plot<-ggplot(data=sub_results, aes(x=Weeks, y=R2, fill=VI)) +
geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7, fill = "lightseagreen") +
theme_bw() +
# theme(legend.position = "bottom", legend.box = "horizontal") +
theme(legend.position = "none") +
guides(colour = guide_legend(nrow = 1)) +
ggtitle(selVI) +
facet_grid(Fieldname~Approach) +
scale_x_continuous(breaks = seq(1, 14, 1)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2, color="black",
position = position_dodge(0), size=2, angle = 90)+
# scale_fill_brewer(palette="Reds") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, face = "bold"),
# axis.title.y = element_blank(),
axis.text = element_text(size = 8, color = "black"),
# axis.text.y = element_blank(),
# axis.ticks.length=unit(-1.5, "mm"),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
bar_plot
# ggsave("VI_performance.pdf", width = 8.5, height = 11, units = "in")
Week 6 and 7 - Three VI + All approaches + All fields
sel_dat <- results_df %>%
filter(Weeks == 6 | Weeks == 7) %>%
filter(VI == "NDVI" | VI == "GNDVI" | VI == "EVI2") %>%
droplevels()
bar_plot<-ggplot(data=sel_dat, aes(x=Weeks, y=R2, fill=VI)) +
geom_bar(stat="identity", position=position_dodge(0.8), width = 0.7) +
theme_bw() +
# theme(legend.position = "bottom", legend.box = "horizontal") +
theme(legend.position = "right") +
guides(colour = guide_legend(nrow = 1)) +
# ggtitle(selVI) +
facet_grid(Fieldname~Approach) +
scale_x_continuous(breaks = seq(1, 14, 1)) +
scale_y_continuous(expand = c(0, 0), breaks = seq(0.0, 1.05, 0.2), limits = c(0.0, 1.05)) +
# geom_text(aes(label=round(R2,2)), vjust=0.5, hjust = -0.2, color="black",
# position = position_dodge(0), size=2, angle = 90)+
# scale_fill_brewer(palette="Reds") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 10, face = "bold"),
# axis.title.y = element_blank(),
axis.text = element_text(size = 8, color = "black"),
# axis.text.y = element_blank(),
# axis.ticks.length=unit(-1.5, "mm"),
plot.title = element_text(size=14, face = "bold")) +
theme(legend.title = element_text(face = "bold", size=8, margin = margin(0,0,7,0)),
legend.text = element_text(vjust = 0.5),
strip.text = element_text(face = "bold", size = 10))
bar_plot
Saving the combined metrics data into CSV